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We examine the validity of the time-dependent Ginzburg-Landau equation of granular fluids for 
a plane shear flow under the Lees-Edwards boundary condition derived from a weakly nonlinear 
analysis through the comparison with the result of discrete element method. We verify quantitative 
agreements in the time evolution of the area fraction and the vertical component to the sheared 
direction of the velocity field, and also found qualitative agreements in the velocity field in the 
sheared direction and the granular temperature. 



I. INTRODUCTION 



Flows of granular particles have been extensively studied due to the importance in powder technology, civil engineer- 
ing, mechanical engineering, geophysics, astrophysics, and applied mathematics and physics 0-Oj]. The characteristic 
properties of the granular flows are mainly caused by the inelastic collisions [5j. In particular, the study of granular 
gases under a plane shear plays an important role in the application of the kinetic theory [6l4l5l|. the shear band in 
a moderate dense flow [ll| [l7( , the long-time tail and the long-range correlations fT§ - [27| , the pattern formation of 
dense flow [28l433| |. the determination of the constitutive equation for dense flow [34 - I36I ] . as well as jamming transition 

The granular hydrodynamic equations based on the kinetic theory well describe the dynamics of moderate dense 
granular gases ja4l5l]. even though its applicability is questionable because of the lack of scale separation and the 
existence of long range correlations, etc. The two-dimensional granular shear flow is an appropriate subject to check 
the validity of the granular hydrodynamic equations, where two shear bands are formed near the boundary and collide 
to form a shear band under a physical boundary condition [I?} ■ A similar shear band is also observed under the 
Lees-Edwards boundary condition. The transient dynamics of the shear band and the hydrodynamic fields can be 
described by the granular hydrodynamic equations, where reasonable agreements with the discrete element method 
(DEM) simulation have been verified [17]. It is also known that a homogeneous state of the two-dimensional granular 
shear flow is intrinsically unstable as predicted by the linear stability analysis [43 - f5(l |. 

To understand the shear band formation after the homogeneous state becomes unstable, we have to develop the 
weakly nonlinear analysis. Recently, Shukla and Alam carried out a weakly nonlinear analysis of the sheared granular 
flow in finite size systems, where they derived the Stuart-Landau type equation for the disturbance amplitude of the 
hydrodynamic fields under a physical boundary condition [5ll - l53l ] . They found the existence of subcritical bifurcation 
in both dilute and dense regimes, while a supercritical bifurcation appears in the medium regime and the extremely 
dilute regime. The Stuart-Landau equation, however, does not include any spatial degrees of freedom and cannot be 
used to study the slow evolution of the spatial structure of shear band. We also notice that the shear rate is fixed to 
unity and cannot be used as a control parameter in their analysis. 

It is also notable that Khain found the coexistence of the solid and liquid phases in his molecular dynamics simulation 
of a dense granular shear flow [HI, [33[ . He demonstrated the existence of a hysteresis loop of the order parameter 
defined as the difference of the densities between the boundary and the center region. It should be noted, however, 
that the mechanism of the subcritical bifurcation based on a set of hydrodynamic equations differs from that observed 
in the jamming transition of frictional particles [43| . 

In our previous work, we have developed the weakly nonlinear analysis for the two-dimensional granular shear flow 
and derived the time dependent Ginzburg-Landau (TDGL) equation for the disturbance amplitude. We introduced a 
hybrid approach to the weakly nonlinear analysis, where the derived TDGL equation is written as a two-dimensional 
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form and has time dependent diffusion coefficients [H^. We have also discussed the bifurcation of the amplitude, 
however, the studies of the numerical solution of the TDGL equation and comparison with the DEM simulation had 
been left as an incomplete part of our previous paper [54] . 

In this paper, we quantitatively examine the validity of the derived TDGL equation for a two dimensional granular 
shear flow from the comparison with the DEM simulation. In Sec. [TTl we review the weakly nonlinear analysis and 
the hybrid approach. In Sec. IIII| we compare the numerical solutions of the TDGL equation with the results of DEM 
simulation. In Sec. IIV1 we discuss and conclude our results. 



II. OVERVIEW OF WEAKLY NONLINEAR ANALYSIS 



In this section, we review our previous results for the weakly nonlinear analysis, where the time evolution for the 
disturbance amplitude is described by the TDGL equation |54[. W e also apply the hybrid approach to the TDGL 
equation to describe the structural changes of the shear band |54| . In Sec. Ill A[ we introduce the basic equations. 
In Sec. Ill B[ we review the weakly nonlinear analysis to derive the TDGL equation. In Sec. Ill CI we derive a 
two-dimensional TDGL equation adopting the hybrid approach to the weakly nonlinear analysis. 



A. Basic equations 



Let us explain our setup and basic equations. To avoid difficulties caused by the physical boundary condition, we 
adopt the Lees-Edwards boundary condition [55| . where the upper and the lower image cells move to the opposite 
directions with a constant speed {7/2. Here, the distance between the upper and the lower image cells is given by 
L. We assume that the granular disks are identical, where the mass, the diameter and the restitution coefficient are 
respectively given by m, d and e. In the following argument, we scale the mass, the length and the time by m, d 
and 2d/U, respectively. Therefore, the shear rate U/L is nondimensionalized as e e 2d/L which becomes a small 
parameter in the hydrodynamic limit L 3> d. 

We employ a set of hydrodynamic equations of granular disks derived by Jenkins and Richman |14j . Although 
their original equations include the angular momentum and the spin temperature, it is known that the spin effects 
are localized near the boundary |56 | and the effect of rotation can be absorbed in the normal restitution coefficient, 
if the friction constant is small [l7l l57l. |58|. Thus, our system is reduced to a system without the spin effects and the 
dimensionless hydrodynamic equations are given by 



(d t + v • V) v = -;A7 . v 
v (d t + v • V) v = -V • P 
0/2) (dt + v- V)0 = -P : Vv 



V-q-x 
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where v % v = (u, w), 9, t and V = (d/d x , d/d y ) are the area fraction, the dimensionless velocity fields, the dimensionless 
granular temperature, the dimensionless time and the dimensionless gradient, respectively. The pressure tensor 
P = (Pij), the heat flux q and the energy dissipation rate x are given in the dimensionless forms as 



P 



q = 

x = 
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respectively, where p(v)6, i(v)9 1 ^ 2 , , q(v)9 1 / 2 , k(u)9 1 / 2 and \(u)9 3 / 2 are the dimensionless forms of the static pressure, 
the bulk viscosity, the shear viscosity, the heat conductivity and the coefficient associated with the gradient of density, 
respectively, and = (VjUj + ViU; — <%V • v)/2 (i,j — x,y) is the deviatoric part of the strain rate tensor. The 
explicit forms of them are listed in Table HI where we adopt the radial distribution function at contact 

I - 7u/16 
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which is only valid for v < 0.7 [59H6 
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p(iz) = ±v[l + (l+e)vg(v)] 

AW = [4(, g (,))^ + 3(l + e)] 



TABLE I: The functions in Eqs.Q-©. 



B. Weakly nonlinear analysis 

To study the slow dynamics of shear band, we need to develop a weakly nonlinear analysis. For this purpose, we 
introduce a long time scale r = e 2 t and long length scales (£, C) = e ( x iV)- We also introduce the neutral solution 
around the most unstable mode q c = (0, q c ) as 

n =A L (C,r)^e^+c.c. , (8) 

where c.c. represents the complex conjugate and </>^_ is the Fourier coefficient at q c . We notice that the amplitude of 
the layering mode A l (^,t) depends only on £ but is independent on £, because any non-layering modes q x j£ are 
linearly stable. Then, we expand A L (£, r) into the series of e as 

A^{C 1 r)=eA\ + e 2 A\ + ^A\ + .. . . (9) 

Substituting Eqs. © and © into the hydrodynamic equations CD)-© an( ^ collecting terms in each order of e we 
obtain an amplitude equation. 

The first non-trivial equation at 0(e 3 ) is the TDGL equation 

d r A\ = a c A\ + Dd\A\ + fiA\\A\\ 2 , (10) 

where D and (3 are listed in Table 2 of Ref. (54j. Here, a c is the maximum growth rate at q c scaled by e 2 . Because of 
the scaling relations D = D and f3 = e/3, we can rewrite the TDGL equation as the equation for the scaled amplitude 
A\ = e x / 2 A\ as 

d T A\ = <r c A\ + DdlA\ + pA\\A\\ 2 . (li) 

It should be noted that the TDGL equation (|10p or (TTTj) can be only used for /3, j3 < 0, i.e., the case of a supercritical 
bifurcation. 

Developing a similar procedure till 0(e 5 ), we also obtain the amplitude equation 

d T A L = a c A L + Dd 2 A L + pA L \A L \ 2 + ejA L \A L \ 4 + 0(e 3 ) , (12) 

where we have introduced A l (C,t) = e^A^r) + eA\{C, T ) + e 2 A\{t,T)] and 7 is also listed in Table 2 of Ref. 
Equation (IT21 can be used for > and 7 < 0, i.e., the case of a subcritical bifurcation. 



C. Hybrid approach to the weakly nonlinear analysis 

Although we derived the TDGL equations (fTTj) and (fT2|) . these equations do not include £ and they are still not 
appropriate to study the two-dimensional structure of shear band. Therefore, we need a new approach, where the 
non-layering mode is coupled with the layering mode. For this purpose, we add a small deviation to the most unstable 
mode as q(r) = q c + <5q(r) and assume 4> n does not change if the deviation <5q(r) is small. Then, Eq. ((5|) can be 
rewritten as 

n ^^ L (C,C,T)^ c e^- z + c.c. , (13) 

where we have introduced z = (£, £) and a ^-dependent amplitude A L (£, £, t). If we also take into account the 
contribution from the non-layering mode, a hybrid solution may be given by 

k = {A L (e,C,T)< c+ ^NL (C5CT)0 NL ) | e , q ( r ). z + c ^ 

c A(e,C,T){# e +< ( L T) }e^)- z + c.c. , (14) 
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where A NL (£, £,t) and 4ffi T ) are the amplitude and the Fourier coefficient of the non-layering mode, respectively. 

Here, we have used a strong assumption that A L (£, T ) an d A NL (£, C t ) are scaled by a common amplitude A(£, £, r) 
in the second line of Eq. (fl4"|) . Expanding j4(£, £ r) as 

A(e,C,r) =eA 1 (C,C,T)+e 2 ^2(e,C ! T)+e 3 A 3 (e,C,r) + ... , (15) 

and carrying out the weakly nonlinear analysis for the hybrid solution <^>h, we found the rescaled amplitude A\ (£, £, r) = 
e 1 / 2 Ai(^, r) for the supercritical bifurcation satisfies 

d T A 1 = a c A x + 5i(r)^iii + D 2 (t)0$0 c A\ + + Ml^il 2 (16) 

at 0(e 3 ), where D\(t) and ^(r) are the time dependent diffusion coefficients. Similarly, we found the higher order 
equation of A(f , £ r) = eV 2 {Ai(f , C, t) + eA 2 (£, (, r) + e 2 A 3 (£, £ r)} as 

r A = a c A + 5i(r)a|i + 5 2 (t)<%<9 c A + 5<9 2 A + ^4|A| 2 + e 7 A|A| 4 + 0(e 3 ) (17) 

for the subcritical bifurcation. The time dependent diffusion coefficients D±(t) and D 2 (t) whose explicit forms are 
given by Eqs. (64) and (65) in Ref. [54[ decay to zero as time goes on. Therefore, Eqs. (fT6)) and (fl7|) are respectively 
reduced to Eqs. (fTT)) and (fT2")l in the long time limit. 

III. DISCRETE ELEMENT METHOD (DEM) SIMULATION 

In this section, we perform the discrete element method (DEM) simulation for a two-dimensional granular shear 
flow to compare the results with the weakly nonlinear analysis. In Sec. IIII A[ we introduce our setup and in Sec. 
IIII B[ we show the time evolution of the density field obtained from the DEM simulation, where the typical transient 
dynamics can be reproduced. In Sec. IIII CI we exhibit the time evolution of the velocity fields and the granular 
temperature, and in Sec. IIII Dl we compare the results of the DEM simulation with the numerical solution of the 
TDGL equation. In the following, we use the same units of mass, length and time as those in the weakly nonlinear 
analysis. 

In Eq. ([15)1. p < for v < 0.245 where the supercritical bifurcation is expected [HI]. If 0.245 < uo_< 0.275, j3 > 
and 7 < 0, thus Eq. (fT7| should be used and the subcritical bifurcation is expected. Unfortunately, /3 > and 7 > 
for vq > 0.275 and neither Eqs. ([Tt?)) nor (|17[) can be used. Therefore, we exhibit our results with i/q = 0.18 for the 
supercritical case, and with vq — 0.26 for the subcritical case, respectively. 

A. Setup 

We adopt the linear spring-dashpot model, where the normal force between the colliding two particles is /„ = 
k n S — r] n S with the overlap S and the relative speed 5. For simplicity, we ignore the tangential contact force, because 
we have already verified the results are unchanged for the realistic value of the friction coefficient by introducing the 
effective restitution coefficient. In our units, the spring and viscosity constants are respectively k n = 500mU 2 /cP 
and rj n — l.OmU/d, and the normal restitution coefficient is e = 0.9. The periodic boundary condition and the 
Lees-Edwards boundary condition with the relative speed U are adopted for the boundaries of the £- and £-axes, 
respectively. Then, we randomly distribute N — 8192 particles in a L* x L* square box with the dimensionless system 
size L* = L/d = 189 (vq = 0.18) and 155 (vq = 0.26), respectively, and randomly distribute the initial velocities 
around the linear velocity profile with the dimensionless shear rate e ~ 10 -2 . 

B. Shear band formation 

Figure[T] (upper panel) displays the time evolution of particles in the DEM simulation for vq = 0.18. Th e hy drody- 
namic fields can be obtained by the coarse graining (CG) procedure developed by Goldhirsch et. al. [63l - t72l |. where 
the CG function is defined as "0( z ) = e _z /ir at z = (£, (). Figure Q] (middle panel) shows the time evolution of the 
area fraction defined as 



N 

^dem(z,t) = ^ s V"(z - z») , 
i=l 



(18) 
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FIG. 1: (Color online) Upper panel: Time evolution of particles in the DEM simulation, where vo = 0.18. Middle panel: Time 
evolution of ^dem(z, t). Lower panel: Numerical solution of Eq. ()16[) . Here, the dimensionless time corresponds to (a) 4.8, (b) 
11.2 and (c) 20.0, respectively. 



where z, = (£j,Ci) an d s = 7r/4 are the dimensionless position of the i-th particle and the dimensionless area of a 
particle, respectively. Figure Q] ( lower panel) shows the numerical solution of Eq. (fT6]) . 

In Fig. [TJ a typical transient dynamics exhibits that (a) the fluctuation with the short wave length is suppressed, 
(b) clusters are generated and merged, and (c) the shear band is generated and the system reaches a steady state. 
Such transient dynamics of shear band is qualitatively similar to the numerical solution of Eq. (|16p. We should 
stress that these results cannot be explained by neither the one-dimensional TDGL equation nor zero-dimensional 
Stuart-Landau equation obtained by the ordinary weakly nonlinear analysis [5ll - [53j . 



C. Velocity fields and granular temperature 



The velocity fields and the granular temperature are defined as 



UDEMIZ, T) = r- , (iyj 

EiW z - z i) 

" DEM(Z ' T) - 2E,^-z t ) ' (20) 
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FIG. 2: (Color online) Upper panel: Time evolution of udem(z, r). Middle panel: Time evolution of u>dem(z, t). Lower panel: 
Time evolution of #dem(z,t). Here, the dimensionless time corresponds to (a) 4.8, (b) 11.2 and (c) 20.0, respectively. 

respectively, where and Vj = Vj — udem ( z i i T ) are the dimensionless velocity of the i-th particle and the di- 
mensionless local velocity, respectively. Figures [2] {upper panel), (middle panel) and (lower panel) display the time 
evolution of udem(z,t), iddem(z,t) and #dem(z,t), respectively, where udem(z,t) and u>dem(z, t) are respectively 
the £ and £ components of udem(z, t). As time goes on, udem(z, t) in the £ direction deviates from the linear profile 
and u>dem(z,t) is almost homogeneous. The time evolution of #dem(z,t) is accompanied with ^dem(z, r), where 
#dem(z,t) is lower in the dense region and higher in the dilute region. 

D. Comparison of the TDGL equation with the DEM simulation 

To test the quantitative validity of the TDGL equation, we compare the numerical solution with the results of DEM 
simulation. At first, we average out z/dem(z,t), «dem(z, t), wdem(z,t) and #dem(z,t) over the £ direction and take 
sample averages from the different 100 time steps. Then, the hydrodynamic fields are written as one-dimensional 
forms ^dem(C: t )) u dem(C,t), w D em(£, t) an d #dem(C,t), respectively. Because ^dem(C t ) an d #dem(£, i~) are 
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approximately symmetric at £ = 0, we introduce 

*dem(C,t) = i{^DEM(C,r)+^DEM(-C,r)} (0 < C < L*/2) , (21) 
DEM (C,r) = i{e DEM (C,T)+0DEM(-C,r)} (0<C<£*/2), (22) 

respectively. On the other hand, the velocity fields are approximately antisymmetric at £ = and we also introduce 

udem(C^) = ^{udem(C,t)-u D em(-C,t)} (0<C<£*/2), (23) 
»dem((,t) = 1;{w D em((,t)~w DE m(-C,t)} (0<(<L*/2), (24) 

respectively. 

In the weakly nonlinear analysis, the hydrodynamic fields are given by the summation of the base state fa — 
(uq, 0, #o) an d the hybrid solution fa. At first, we project fa on the £-axis as 

^(C,T)~I(C,T)^ e ^Mf + c.c. , (25) 

where qc{r) = q c — r is the C component of q(r) [73[ and we ignore because </>q(^) exponentially decays to 

zero and the following results are unchanged even if we take into account ^q^)- We note that tp^ c is defined as 
<fiq c = (i , qc ,iu qc ,iw qc ,9 qc ) T with the imaginary unit i, where v Qc , u qc , w qc and 9 Qc have been given in our previous 
study [54j |. If we ignore the higher order terms in Eq. (fl~5|) . A((, r) may be given by the numerical solution of Eq. 
(ITBI projected on the £-axis. Then, the hydrodynamic fields are given by <^tdgl(C t ) — fa + 01i(C t )j where each 
component of 0tdgl(C t ) is written as 

zTdgl(C,t) = ^ + 2^ c I(C,T)cos(g c (r)C) , (26) 

utdgl(C.t) = C-2u 9c I(C,T)sin(< Zc (r)C) , (27) 

«Mt>gl(C.t) = -2^ c I(C,T)sin(g c (r)C) , (28) 

0tdgl(C,t) = 9 Q + 20 q J(C,T)cos(q c (T)C) , (29) 

respectively, where the factor 2 comes from the complex conjugate. 

Figures [3] and S] display the time evolution of the hydrodynamic fields for the supercritical case (vo = 0.18) and the 
subcritical case (vo = 0.26), respectively, where the symbols represent Eqs. (l2~Tj) - (l24l obtained by the DEM simulation 
and the lines represent the scaling functions 

^TDGL(C,T) = a^XTDGL(C/Ci(T),T/T*) (X = V, U, W,9) (30) 

with the scaling factors a* x , Cx( T ) an( i T *i respectively. In these figures, S > tdgl(Ci t ) an d ^tdgl(Cj t ) are quantita- 
tively agreed with Pdem(Cj t ) an( i ^dem(C> t )j respectively. We should note that we could not get any reasonable 
agreement between the £ component of the velocity field in a numerical solution of full set of the granular hydro- 
dynamic equations and the result of DEM simulation in our previous work [l7| . We can also see the qualitative 
agreements in the £ component of the velocity field in the the granular temperature. 



IV. DISCUSSION AND CONCLUSION 



In this paper, we examine the validity of the TDGL equation for a two-dimensional sheared granular flow from 
the comparison with the results of the DEM simulation with the aid of the CG method. The results of the TDGL 
equation, at least, qualitatively agree with the results of the DEM simulation. Such transient dynamics cannot be 
reproduced by neither the one dimensional TDGL equation nor the zero dimensional Stuart-Landau equation derived 
by the ordinary weakly nonlinear analysis. We also obtain that the velocity fields and the granular temperature 
qualitatively agree with the solution of the TDGL equation. 

We compare the one dimensional hydrodynamic fields obtained by averaging the results of the DEM simulation 
with the scaled forms of the numerical solution of the TDGL equation, where we find the qualitative agreements 
in the area fraction and the velocity field in the ( direction, while there are still deviations in the velocity field in 
the £ direction and the granular temperature. In our previous work, we demonstrated that the hydrodynamic fields 
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FIG^3: Time evolution of (a) voem{C,,t) and ^tdgl(C,''~), (b) ^DEM(C r ) and mtdglCC, 1 "), ( c ) w em{C,,t) and wtdglCC 1 "), 
(d) #DEM(Ci r ) an( i $tdgl(C,t), respectively, for the supercritical case (fo = 0.18), where the solid squares and the solid lines 
correspond to the dimensionless time 4.8, the solid circles and the hashed lines correspond to the dimensionless time 11.2, and 
the solid triangles and the dotted lines correspond to the dimensionless time 20.0, respectively. Here, we have introduced the 
scaling factors r* ~ 0.14, a* ~ 0.24, a* ~ 0.02, aj, ~ 1.96 and ag ~ 0.02, respectively, and we also use Ct( T ) — 1-6,0.9,0.75, 
Cu(t) ~ 0.8,0.8,0.8, Cw(t) ~ 1.5, 1.1, 1.8 at the the dimensionless time 4.8, 11.2 and 20.0, respectively, and Ce( r ) — 1-6, 1-35 
at the dimensionless time 11.2 and 20.0, respectively. It should be noted that we do not show the result of the granular 
temperature at the dimensionless time 4.8, because it homogeneously distributed around 9o and the fluctuation is too large to 
plot in the same figure. 



obtained by the DEM simulation can be reasonably explained by the numerical solution of the full hydrodynamic 
equations by Jenkins and Richmann except for w(z,t) [lil . Il5l . [l7| . In the present work, even though we need to 
introduce the scaling factors, the results of the DEM simulation is qualitatively reproduced by the numerical solution 
of the TDGL equation. It is needless to say that more precise analyses to remove the scaling factors will be important 
and is left to our future works. 

In conclusion, the numerical solution of the TDGL equation can qualitatively explain the time evolution of the 
hydrodynamic fields obtained by the DEM simulation. 
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